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Abstract. Linear fluctuating hydrodynamics is a useful and versatile tool for describing 
fluids, as well as other systems with conserved fields, on a mesoscopic scale. In one 
spatial dimension, however, transport is anomalous, which requires to develop a nonlinear 
extension of fluctuating hydrodynamics. The relevant nonlinearity turns out to be the 
quadratic part of the Euler currents when expanding relative to a uniform background. 
We outline the theory and compare with recent molecular dynamics simulations. 
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1 Introduction, long time tails for simple fluids 

In the mid 1950ies, Green [1] and Kubo [2, 3] discovered that transport coefficients for 
simple fluids can be obtained through a time-integral over the respective total current 
correlation function. For tracer diffusion such a connection is more immediate and was 
understood much earlier. But the then novel insight was that collective transport co¬ 
efficients, such as viscosity and thermal conductivity, follow the same pattern. Thus it 
became a central issue to determine the time decay of such current correlations. With 
essentially no tools available this amounted to an impossible task. The static equilibrium 
correlations were known to decay exponentially fast, as conhrmed by a convergent series 
expansion. But for the dynamics one would have to deal with a huge set of coupled 
differential equations. At the time the only theoretical tool available was the Boltzmann 
equation valid at low density. Kinetic theory predicts an exponential decay for the current 
time correlations and it was tacitly assumed that such behavior would extend to moderate 
densities. Alder and Wainwright [4] tried to check the situation in a pioneering molecular 
dynamics simulation of 500 hard disks, resp. hard spheres, at periodized volume 2, 3, 
and 5 times larger than close packing. For tracer diffusion they convincingly observed a 
power law decay as dimension d = 2, 3, which was baptized “long time tail”. They 

also argued that the same behavior should hold for collective transport. Theory quickly 
jumped in and predicted a decay as for viscosity and thermal conductivity [5, 6]. 

There are several theoretical schemes and they all arrive at the same prediction, which 
of course increases their conhdence level. As recognized for some time, the most direct 
approach is linear fluctuating hydrodynamics plus small nonlinear perturbations. We 
refer to the recent monograph [7] for a comprehensive discussion. Here I provide only a 
rough sketch of the method with the purpose to explain why one dimension is so special. 
In the physical dimension a simple fluid has hve conservation laws and correspondingly 
fluctuating hydrodynamics has to deal with a hve component held, where the momentum 
components are odd, density and energy are even under time reversal. As well known 
[8, 7], the full structure has to be used in order to arrive at quantitative predictions. But 
the argument becomes even more direct for the (unphysical) case of a single conservation 
law. 

Let us thus consider the scalar held p{x,t), which for concreteness is called density. 
Space is a; G and time is t G M. p{x, t) is a huctuating held. On the macroscopic scale 
huctuations are not visible and p satishes the hyperbolic conservation law 

+V ■ J(p) = 0, (1.1) 

where j is the density current. To also include mesoscopic details, in particular to in¬ 
corporate huctuations, one argues that the current has, in addition to the deterministic 
part, also a random contribution which is essentially uncorrelated in space-time. Since 
huctuations are always associated with dissipation, on a more rehned scale Eq. (1.1) 
becomes 

dtp + ^ ■ DWp + ai) =Q. ( 1 . 2 ) 
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( 1 . 3 ) 


Here ^{x,t) is Gaussian white noise with mean zero and covariance 


= 6aa'S{x - x')6{t - t') , tt, «' = 1, 


a is the noise strength and D is the diffusion constant. They are both treated as numbers. 
Physically they will depend on the density. But this would be higher order effects, which 
are ignored in our discussion. 

The goal is to compute the density time correlations in the stationary regime, which 
is no easy task, since (1.2) is nonlinear. But correlations can be thought of as imposing 
at t = 0 a small density perturbation at the origin and then record how the perturbation 
propagates in space-time. For this purpose one might hope to get away with linearizing 
(1.2) at the uniform background density po as po + g{x,t), which yields 

+ DVg + ai) =0, (1.4) 


where ' refers to differentiation w.r.t p. (1.4) is a linear Langevin equation, hence solved 
easily. Since p is a deviation from the uniform background, we are interested in the space- 
time stationary process with zero mean. First note that (1.4) has a unique time-stationary 
zero mean measure, which is Gaussian white noise in the spatial variable, 

{q{x)q{x')) = /2D)6{x — x'). (1.5) 

For the stationary space-time covariance one obtains 

{q{x, t)p(0, 0)) = {o^l2D)p{x - ct, t) , (1.6) 

where p is the Gaussian transition kernel, 

p{x, t) = ( 47 rZl|t |)-''/2 ^ _ xyAD\t\) , (1.7) 

and c = j'{po)- By (1.4) the fluctuating current is given by 

J = eg — DVg + a^. (1.8) 


For the stationary correlations of the total current one arrives at 


dx{J{x, t) ■ 57(0, 0)) = da‘^6{t). 


( 1 . 9 ) 


No surprise. A density fluctuation propagates with velocity c and spreads diffusively. 
The currents are delta-correlated, which should translate into exponential decay for the 
underlying microscopic system. But before jumping at such conclusions one has to study 
the stability of (1.9) against including higher orders in the expansion. By power counting 
the next to leading term is the nonlinear current at second order, which amounts to 

dtg V ■ (cp -|- Gg^ — Z)V g ) = 0 , G = ' (po) • (1-10) 

The task is to compute the current correlation (1.9) perturbatively in G. 
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We first remove c by switching to a moving frame of reference. Secondly we note 
that, quite surprisingly, white noise is still time-stationary under the nonlinear Langevin 
equation (1.10). Since, as argued above, spatial white noise is time-stationary under 
the linear part, one only has to check its invariance under the deterministic evolution 
dtQ = —G ■ Vp^. Formally its vector held is divergence free, since one can choose a 
divergence free discretization [9], compare with the discussion in Section 3 below. Hence 
the “volume measure” is preserved and one has to check only the time change of the 
logarithm of the stationary density, 

— [ dx^g{xy = — f dxQ{x)G ■ Wg‘^{x) = ^ f dxG ■ Wg^{x) = 0 . (l-H) 

dt jRd jRd 

For an expansion in G ■ it is most convenient to use the Fokker-Planck generator, 
denoted by L = Lq + Li, where Lq corresponds to the Gaussian part and Li to the 
nonlinear how. We dehne 


t) = {g{x, t)g( 0 , 0)) = {e(x)e‘-‘g( 0 )) , 


( 1 . 12 ) 


average with respect to the time-stationary Gaussian measure, and plan to use the general 
second moment sum rule 


^ f dxx^S{x,t) = [ dx{J{x,t) ■ J{0,0)), 

dt Jmd ./TOd 


( 1 . 13 ) 


which follows directly from the conservation law. By second order time-dependent per¬ 
turbation theory, 

S{x,t) = {p(a;)e^°*p(0)) + f dti(p(a;)e^°(‘-*^)Lie^“‘ip(0)) 

Jo 

+ f dt 2 r + 0{\Gf). (1.14) 

Jo Jo 

Log is linear, while Lip is quadratic in p. Thus the second term on the right is odd in p 
and vanishes. For any functional, F, 


{g{x)LiF{g)) = -{{Lig{x))F{g )), 


(1.15) 


see (4.8), and the left Li is swapped over to act on q{x). By translation invariance, 

G ■ V can be pulled in front as acting on x. Gombining the terms one arrives at 

S{x,t) = {g{x)e^°^g{0)) + f dt 2 f dti(G-V)^ 

Jo Jo 

X I dxi I dx2p{x - X2,t - t2)p{xi,ti){g{x2)‘^e^°d2-ti) (1-16) 


The Gaussian average is computed as 

{ 0 {x 2 fe^°^g{xif) = 2 {{a^/ 2 D)p{x 2 - Xi, t)Y + s-c . 


(1.17) 
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The self-contraction does not depend on x, hence vanishes when applying G ■ V. We 
insert in (1.13). Working out the integrals yields, including second order, 


dx(J(x, t) ■ J(0, 0)) = dcr''‘6(t) + 4(x^/2D)^lGl^p(0,2t), 


( 1 . 18 ) 


which is the claimed long time tail for a scalar conserved held in d dimensions. 

Eq. (1.10) is singular at short distances, the worse the higher the dimension. In 
physical systems there is a natural cut-off at the microscopic scale which is simply ignored 
in (1.10). There are many possibilities to improve, but the serious constraint is to obtain 
a still manageable nonlinear stochastic equation. The noise should remain 5-correlated 
in time so to preserve the Markov property of the time evolution. One could smoothen 
in space, but thereby may loose the information on the time-stationary measure. To my 
experience the best compromise is to adopt the obvious spatial discretization. Then one 
has a set of coupled stochastic differential equations. For them one can check rigorously 
the time-stationary measure and identities as e.g. (1.15). On the perturbative level 
this then leads to the continuum expressions as (1.16). For a simple fluid the current 
correlations are bounded and will not diverge near t = 0. So in (1.18) only the long time 
prediction can be trusted. 

What can be learned from the long time tails? In dimension d > 3 the decay is 
integrable. Thus, as assumed implicitly beforehand, the model has a hnite diffusivity. 
Higher order terms in the expansion will modify prefactors but should not alter the 
exponent for the time decay. Dimension d = 2 is marginal. The diffusivity is weakly 
divergent. In principle, say, a system of hard disks has inhnite viscosity and thermal 
conductivity. But since the divergence is only logarithmic one could convert the result 
into a weak system size dependence. In one dimension the conductivity is truly inhnite. 
Obviously, not even the power law as based on the perturbative expansion (1.16) can be 
trusted and one needs to develop non-perturbative techniques. 

Eq. (1.10) for d = 1 is known as stochastic Burgers equation, which we record for 
later reference as 

dtQ — dx{\XQ^ +'odxQ + Q , xeM. (1-19) 

One can introduce a potential through q = dxh. Then h satishes the one-dimensional 
version of the Kadar-Parisi-Zhang (KPZ) equation [10], 

dth = ^X{dxhy + udlh + , ( 1 - 20 ) 

for the moment using the more conventional symbols for the coefficients. Over the last 
hfteen years many properties, including rigorous results, of the KPZ equation have been 
obtained. While this is not the place to dive into details, I note that the solution is 
continuous in x, but so singular that (d^h)^ is ill-dehned. However it is proved that the 
ultraviolet divergence is very mild and can be tamed by what would be an inhnite energy 
renormalization in a 1-1-1 dimensional quantum held theory, compare with the discussion 
in [11]. More precisely one chooses a regularized version of (1.20) by replacing ^{x,t) 
through the spatially smoothed version t) = f (p(x — x')^(x', t)dx' with (p > 0, even, 
smooth, of rapid decay at inhnity, and normalized to 1. It can be proved that then the 
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solution to (1.20) is well-defined. One introduces an ultraviolet cut-off of spatial size e by 
choosing the ^-function sequence = e“^(p(e“^x) and substituting the noise ^ by its 

smoothed version Let us denote the corresponding solution of the KPZ equation 

by he(x, t). Then h^(x, t) —vj, i.e. he{x, t) viewed in the frame moving with the diverging 
velocity f (p(x)^dx, has a non-degenerate limit as e —)■ 0, which coincides with 

the Cole-Hopf solution of the KPZ equation [12, 13]. 

In this context Hairer recently developed a solution theory, for which he was awarded 
the 2014 Fields Medal. His theory works for the KPZ equation, as well as a large class 
of other singular stochastic partial differential equations, and for general cutoffs [14, 15]. 
Of course, the solution theory studies the small scale structure of solutions, and not the 
large scale, where the universal behaviors of interest are observed. 

There is an exact formula for S{x,t), which involves Fredholm determinants. In the 
long time limit 

S{x,t) ~ {a'^/2u){V2X\t\)~^^^fKPz{V2X\t\)~^^^x) . (1.21) 

The scaling function /kpz will reappear below, where more details are given. The Burgers 
current reads 

J{x,t) =-{^Xg^ + uda^g + a^) (1.22) 

and for its total correlation function one obtains, using the sum rule (1.13), 

[ dx{J{x,t)J{0,0)) ({a‘^/2iy){y/2xy [ dx^x"^fKPz{x)){V2Xt)~‘^^^ (1.23) 

valid for large t. Note that the true decay turns out to be to be contrasted with 

the perturbative result In fact the power 2/3 was argued already in [5] using a self- 

consistent scheme, see also [16]. As a fairly unusual feature, the non-universal coefficients 
are given directly in terms of the parameters of the stochastic Burgers equation. In fact, 
the particular form can be guessed from the scaling properties of Eq. (1.19). Only for 
the scaling function /kpz and the pure number a/2 one has to rely on the exact solution, 
which is the result of an intricate analysis using lattice type models [17, 18, 19, 20], 
replica computations for the KPZ equation [21, 22], and the finite time exact solution of 
the stationary KPZ equation itself [23]. 

In these notes we will explain how nonlinear fluctuating hydrodynamics, in the same 
spirit as already explained for a scalar field in one dimension, can be used to predict the 
asymptotic form of the equilibrium space-time correlations of anharmonic chains. Chains 
are one-dimensional objects and it may seem as if we have accomplished the task already. 
Well, we have not even written down a Hamiltonian. So first we have to dwell on a few 
general properties of anharmonic chains. In particular we will see that they have three 
conserved helds, generically. The corresponding Euler equations are derived, thereby iden¬ 
tifying the macroscopic currents. But now we are forced to handle several conserved helds. 
While it is not so difficult to write down the multi-component generalization of (1.19), 
the analysis of the Langevin equation will be more complicated with no exact solution at 
help. A second major task will be to test the quality of these predictions by comparing 
with molecular dynamics simulations. 
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2 Anharmonic chains 

Conservation laws and equilibrium time correlations. We consider a classical flnid 
consisting of particles with positions qj and momenta pj, j = 1,..., N, qj,Pj G M, possible 
bonndary conditions to be discnssed later on. We nse nnits snch that the mass of the 
particles eqnals 1. Then the Hamiltonian is of the standard form, 

N N 

< = E k? + i E • (2-1) 

j=l i¥=j=^ 

with pair potential V{x) = V{—x). The potential may have a hard core and otherwise 
is assnmed to be short ranged. The dynamics for long range potentials is of independent 
interest [24] , bnt not discnssed here. Fnrthermore the potential is assnmed to be thermo¬ 
dynamically stable, meaning the validity of a bonnd as ~ Qj) ^ ^ ~ BN for 

some constants A and R > 0. A snbstantial simplihcation is achieved by assnming a hard 
core of diameter a, i.e. R(a:) = oo for |a;| < a, and restricting the range of the smooth 
part of the potential to at most 2a. Then the particles maintain their order, qj < qj+i, 
and in addition only nearest neighbor particles interact. Hence Hff simplihes to 

N N-1 

= E id + E - ''j) ■ (2,2) 

j=l j=l 

As a, at hrst sight very different, physical realization, we conld interpret if at as describing 
particles in one dimension conpled throngh anharmonic springs which is then nsnally 
referred to as anharmonic chain. 

In the second interpretation the spring potential can be more general than anticipated 
so far. No ordering constraint is reqnired and the potential does not have to be even. 
To have well dehned thermodynamics the chain is pinned at both ends as gi = 0 and 
qN+i = (-N. It is convenient to introdnce the stretch Vj = — qj. Then the bonndary 

condition corresponds to the microcanonical constraint 

N 

E ■ 

i=i 
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(2.3) 


Switching to canonical equilibrium according to the standard rules, one the arrives at the 
obvious condition of a hnite partition function 

(2.4) 

using the standard convention that the integral is over the entire real line. Here /3 > 0 is 
the inverse temperature and P is the thermodynamically conjugate variable to the stretch. 
By partial integration 


P = -Z{P, /3)-^ J dxV'ix) , (2.5) 

implying that P is the average force in the spring between two adjacent particles, hence 
identihed as thermodynamic pressure. To have a hnite partition function, a natural 
condition on the potential is to be bounded from below and to have a one-sided linear 
bound as V{x) > oq -l- for either a; > 0 or x < 0 and bo > 0. Then there is a 
non-empty interval I{(3) such that Z{P,f3) < cx) for P G /(/3). For the particular case of 
a hard-core huid one imposes P > 0. 

Note: The sign of P is chosen such that for a gas of hard-point particles one has the 
familiar ideal gas law P = l/(3i. The chain tension is —P. 

Famous examples are the harmonic chain, Vha(a^) = the Fermi-Pasta-Ulam (FPU) 
chain, Vfpu(3^) = -|- \(3x*, in the historical notation [25], and the Toda chain 

[26], V{x) = e“^, in which case P > 0 is required. The harmonic chain, the Toda chain, 
and the hard-core potential, Vhc(x) = oc for |x| < a and Vhc{x) = 0 for |x| > a, are in 
fact integrable systems which have a very different correlation structure and will not be 
discussed here. Except for the harmonic chain, one simple way to break integrability is 
to assume alternating masses, say rrij = mo for even j and ruj = rui for odd j. 

We will mostly deal with anharmonic chains described by the Hamiltonian (2.2), in¬ 
cluding one-dimensional hard-core fluids with a sufficiently small potential range. There 
are several good reasons. Firstly on the level of fluctuating hydrodynamics a generic 
one-dimensional fluid cannot be distinguished from an anharmonic chain. Thus with 
the proper translation of the various terms we would also predict the large scale corre¬ 
lation structure of one-dimensional fluids. The second reason is that in the large body 
of molecular dynamics simulations there is not a single one which deals with an “hon¬ 
est” one-dimensional fluid. To be able to reach large system sizes all simulations are 
performed for anharmonic chains. In addition, from a theoretical perspective, the equilib¬ 
rium measures of anharmonic chains are particularly simple in being of product form in 
momentum and stretch variables. Thus material parameters, as compressibility and sound 
speed, can be expressed in terms of one-dimensional integrals involving the Boltzmann 
factor q-PN{^)+p^) ^ U(x), and x. 

Anharmonic chains should be thought of as a particular class of 1-1-1 dimensional 
held theories. Thus qj is viewed as the displacement variable at lattice site j and not 
necessarily the physical position of the j-th particle on the real line. There is a simple 
translation between both pictures, but we will stick to the held theory point of view. For 
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a fluid with unlabeled particles and, say, a bounded potential, the equivalence is lost and 
only the fluid picture can be used. 

This being said, we follow the standard rules. We write down the dynamics of stretches 
and momenta and identify the conserved helds. From there we infer the microcanonical 
and canonical equilibrium measures. For slowly varying equilibrium parameters we de¬ 
duce the Euler equations. In particular, their version linearized at uniform equilibrium 
will constitute the backbone in understanding the equilibrium time correlations of the 
conserved fields. 

The dynamics of the anharmonic chain is governed by 

= Pi ) j^Pj = V'iQj+i - Qj) - V'iQj - Qj-i) ■ ( 2 . 6 ) 

For the initial conditions we choose a lattice cell of length N and require 

Qj+N = Qj + iN , Pj+N =Pj (2.7) 

for all j G Z. This property is preserved under the dynamics and thus properly mimics a 
system of hnite length N. The stretches are then iV-periodic, tj+n = rj, and the single 
cell dynamics is given by 

= pj+i - pj , ^pj = V\rj) - , ( 2 . 8 ) 

j = 1,... ,N, together with the periodic boundary conditions pi+N = Pi, and 

the constraint (2.3). Through the stretch there is a coupling to the right neighbor and 
through the momentum a coupling to the left neighbor. The potential is dehned only up 
to translations, since the dynamics does not change under a simultaneous shift of V{x) 
to V{x — a) and Vj to Vj + a, in other words, the potential can be shifted by shifting 
the initial r-£eld. Note that our periodic boundary conditions are not identical to fluid 
particles moving on a ring, but they may become so for large system size when length 
fluctuations become negligible. Both equations are already of conservation type and we 
conclude that 


d A d A 

i=i j=i 

(2.9) 

We define the local energy by 

^i = \P] + ^ (^i) • 

(2.10) 

Then its local conservation law reads 


= Pj+iV'{rj) - PjV'{rj_i) , 

(2.11) 

implying that 

cl \^ 

1=1 

(2.12) 
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At this point we assume that there are no further local conservation laws. Unfortu¬ 
nately our assumption, while reasonable, is extremely difficult to check. It certainly rules 
out the integrable chains, which have N conservation laws. There is no natural example 
known with, say, seven conservation laws. But the mere fact that there are exceptions 
implies that close to integrability the predictions from fluctuating hydrodynamics could 
be on time scales which are not accessible. The parameters entering in fluctuating hy¬ 
drodynamics depend smoothly on the potential. Thus as one approaches, for example, 
the Toda potential no abrupt changes will be detected. In this sense, the theory cannot 
distinguish the Toda chain from a FPU chain both at moderate temperatures. There is 
another limitation which can be addressed more quantitatively. If V (x) -|- Px has a unique 
minimum, then at very low temperatures the potential is close to a harmonic one. This 
feature will be properly reflected by fluctuating hydrodynamics through the temperature 
dependence of the coupling coefficients and the noise strength. 

The microcanonical equilibrium state is dehned by the Lebesgue measure constrained 
to a particular value of the conserved helds as 

N N N 

Y^r, = tN, ^ (Ip] + V(rj)) = cN (2.13) 

j = l j = l j = l 

with i the stretch, u the momentum, and e the total energy per particle. In our context 
the equivalence of ensembles holds and computationally it is of advantage to switch to the 
canonical ensemble with respect to all three constraints. Then the dual variable for the 
stretch £ is the pressure P, for the momentum the average momentum, again denoted by 
u, and for the total energy e the inverse temperature (3. For the limit of inhnite volume 
the symmetric choice j G [—A^, ...,iV] is more convenient. In the limit iV —)• cx) either 
under the canonical equilibrium state, trivially, or under the microcanonical ensemble, by 
the equivalence of ensembles, the collection {rj,pj)j^z are independent random variables. 
Their single site probability density is given by 

_ ^2.14) 

Averages with respect to (2.14) are denoted by (■)p,/ 3 ,u- The dependence on the average 
momentum can be removed by a Galilei transformation. Hence we mostly work with 
u = 0, in which case we merely drop the index u. We also introduce the internal energy, 
e, through e = -|- e, which agrees with the total energy at u = 0. The canonical free 

energy, at u = 0, is dehned by 

GiP,f3) = -/3-i( - ilog/3 + logZ(P,/3)) . (2.15) 

Then 

£ = {ro)p ,0 , e = d^{(3GiP, (3)) - P£ = ^ + {V{ro))p,p . (2.16) 

The relation (2.16) dehnes {P, (3) h-)■ {£{P,/3), e{P, (3)), thereby the inverse map (-^,6) h->■ 

{P{£,e), /3{£,e)), and thus accomplishes the switch between the microcanonical thermo¬ 
dynamic variables £, e and the canonical thermodynamic variables P, (3. 
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It is convenient to collect the conserved fields as the 3-vector g = {gi,g 2 ,g 3 ), 

9i.ht) = , (2.17) 

0) = diJ)- Then the conservation laws are combined as 

+ J{j + l,t) - J{j,t) = 0 (2.18) 

with the local current functions 

JU) = ( - P„ -V(o-i). -P,y'l.r,-i)). (2.19) 

Our prime interest are the equilibrium time correlations of the conserved fields, which are 
defined by 


Saa'{j,t) — {9a{j,t)9a'i0,0))p^p — (S'o (0, 0))p^_a(5fQ,'(0, 0))p_/3 , (2.20) 

a, a' = 1, 2, 3. The infinite volume limit has been taken already and the average is with 
respect to thermal equilibrium at u = 0. It is known that such a limit exists [28]. Also the 
decay in j is exponentially fast, but with a correlation length increasing in time. Often it 
is convenient to regard S{j, t), no indices, as a 3 x 3 matrix. In general, S{j, t) has certain 
symmetries, the first set resulting from space-time stationarity and the second set from 
time reversal, even for a = 1,3, odd for a = 2, 

5'«a'(j) t) = Sa'a{-j, -t) , S^a'U, t) = (-1)“+“'S'„„'(j, -t) . (2.21) 

At f = 0 the average (2.20) reduces to a static average, which is easily computed. The 
fields are uncorrelated in j, i.e. 

S{j,0)=5joC (2.22) 

with the static susceptibility matrix 

/(^o;^o)p,/3 0 (’"o;To)p,/3 \ 

0=0 0 . (2.23) 

\(’"o; To)p,/ 3 0 + (fo; fo)p,/3/ 

Here, for X,Y arbitrary random variables, {X]Y) = (XY) — {X){Y) denotes the second 
cumulant and Vq = V{ro), following the same notational convention as for cq. The 
conservation law implies the zeroth moment sum rule 

(j, = (t 0) = • (2.24) 

j&Z j&Z 

An explicit computation of S{j,t) is utterly out of reach. But with the current com¬ 
puter power MD simulations have become an essential source of information. A broader 
coverage will be provided in Section 5. Just to have a first impression I show in Fig. 1 
a recent MD simulation of the correlator for a FPU chain. One notes the central peak. 
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called heat peak, which is standing still, and two symmetrically located peaks, called 
sonnd peaks which move ontwards with the speed of sonnd. The peaks broaden with a 
certain power law which will have to be discnssed. One expects, better hopes for, self¬ 
similar shape fnnctions, at least for snfhciently long times. We do not know yet their 
form. Bnt the central peak seems to have fat tails while the sonnd peaks fall off more 
rapidly, at least towards the ontside of the sonnd cone. The area nnder each peak is 
preserved in time and normalized to 1 in onr plot. If the chain is initially pertnrbed near 
0 and the response in one of the conserved fields is observed at (j, t), then one records a 
signal, which is a linear combination of the peaks in Fig. 1, the coefficients depending on 
the initial pertnrbation. It might happen that one peak is missing. If the pertnrbation 
is orthogonal to all three physical fields, then there is no peak at all, only low amplitude 
noise. 



Figure 1: Heat peak and sound peaks, area normalized to 1, at times t = 800,1300, 2700, 
for a FPU chain with potential parameters a = 2, /3 = 1, pressure P = 1, and temperature 
/3-1 = 0.5. 

Our goal is to predict the large scale structure of the correlator Saa' (j, t) encoding the 
propagation of local perturbations of the equilibrium state. On the crudest level, they 
should be captured by linearized hydrodynamics, to which we turn next. 

Linearized hydrodynamics. We start the dynamics from a product measure of the form 
(2.14), but replace the uniform P, u,(3 by slowly varying spatial helds P{ej), u{ej), (3{ej) 
with e -C 1. e~^ is the macroscopic scale measured in lattice units. Equivalently we may 
regard e as the lattice spacing. The relation (2.16) then dehnes also the slowly varying 
helds u(ej), e(ej). Because of the conservation laws the time change of such a state is 
slow and varies only over microscopic times of order e~^t with macroscopic t of order 1. We 
average Eq. (2.18) over the slowly varying initial state, which gives then the time change 
of the average locally conserved helds. The expectation of the currents is more difficult. 
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The difference in j becomes ed^ on the macroscopic scale. Since the currents are functions 
depending only on at most two neighboring lattice sites, to lowest approximation their 
average can be computed in the equilibrium state with the corresponding local values of 
the fields, i.e. £(ej, u(ej, c(ej, Therefore we define the hydrodynamic 

Euler currents by the equilibrium averages 

(i(j))qu,e = ( - U, P(£, e - |u2), UP(£, e - ^u^)) = u, e) (2.25) 

with P{i, e) defined implicitly through (2.16). Our argument then leads to the macroscopic 
Euler equations 


dti — dxU = 0 , dfU + dxP{i,t — ^u'^) = 0 , dtt + dx:[uP{i,c — = 0 . (2.26) 

We refer to a forthcoming monograph [28], where the validity of the Euler equations is 
proved up to the first shock. Since, as emphasized already, it is difficult to deal with 
deterministic chaos, the authors add random velocity exchanges between neighboring 
particles which ensure that the dynamics locally enforces the microcanonical state. 

We are interested here only in small deviations from equilibrium and therefore linearize 
the Euler equations as i + Ui{x), 0 + U 2 {x), e + U 3 (x) to linear order in the deviations 
u(x). This leads to the linear equation 

dtu{x,t) + dxAu{x,t) = 0 (2.27) 

with 

/ 0 -1 0 \ 

A = \diP 0 ^eP . (2.28) 

\ 0 P 0 / 

Here, and in the following, the dependence of H, C and similar quantities on the back¬ 
ground values £, u = 0, e, hence on P,/3, is suppressed from the notation. Beyond (2.24) 
there is the first moment sum rule which states that 

J2jS{jA)=ACt. (2.29) 

iez 

A proof, which in essence uses only the conservation laws and space-time stationarity of 
the correlations, is given in [29], see also see [30, 31]. Microscopic properties enter only 
minimally. However, since C = and Eq. (2.29) implies the 

important relation 

AC = {ACf = CA'^ , (2.30) 

with denoting transpose. Of course, (2.30) can be checked also directly from the defini¬ 
tions. Since O > 0, A is guaranteed to have real eigenvalues and a nondegenerate system 
of right and left eigenvectors. For A one obtains the three eigenvalues 0, ±c with 

c2 = -dtP + Pd^P > 0 . (2.31) 
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Thus the solution to the linearized equation has three modes, one standing still, one right 
moving with velocity c and one left moving with velocity —c. Hence we have identihed 
the adiabatic sound speed as being equal to c. 

(2.27) is a deterministic equation. But the initial data are random such that within 
our approximation 

{Ua{x, 0)Ua'{x\ 0)) = Cq,q/(5(x — x'). (2.32) 

To determine the correlator S{x,t) with such initial conditions is most easily achieved by 
introducing the linear transformation R satisfying 

RAR~^ = diag(—c, 0, c), RCR^ = 1. (2.33) 


Up to trivial phase factors, R is uniquely determined by these conditions. Explicit for¬ 
mulas are found in [29]. Setting cj) = Au, one concludes 

d(j)a + Cadx(j)a = 0 , Q! = -1, 0, 1 , (2.34) 

with c = (—c, 0,c). By construction, the random initial data have the correlator 

{(j)a{x, 0)(j)a'{x', 0)) = 6aa'5{x - x') . (2.35) 

Hence 

{(j)a{xR)(j)a'{0,0)) = 6aa'S{x - Cat) . (2.36) 

We transform back to the physical helds. Then in the continuum approximation, at the 
linearized level, 

S{x, t) = i?“^diag(5(x -|- ct), S{x), 6{x — ct))R~^ (2.37) 

with R-^ = (R-^)^. 

Rather easily we have gained a crucial insight. S{j,t) has three peaks which separate 
linearly in time. For example, t) has three sharp peaks moving with velocities ±c, 0. 
Physically, one expects such peaks not to be strictly sharp, but to broaden in the course 
of time because of dissipation. This issue will have to be explored in great detail. It 
follows from the zeroth moment sum rule that the area under each peak is preserved in 
time and thus determined through (2.37). Hence the weights can be computed from the 
matrix R~^, usually called Landau-Plazcek ratios. A Landau-Placzek ratio could vanish, 
either accidentally or by a particular symmetry. An example is the momentum correlation 
S 22 {j,t)- Since (R“^) 2 o = 0 always, its central peak is absent. 

For integrable chains each conservation law generates a peak. Thus, e.g., Sii{j,t) of 
the Toda chain is expected to have a broad spectrum expanding ballistically, rather than 
consisting of three sharp peaks. 


3 Nonlinear fluctuating hydrodynamics 

Euler currents to second order. The broadening of the peaks results from random 
fluctuations in the currents, which tend to be uncorrelated in space-time. Therefore the 
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crudest model would be to assume that the current statistics is space-time Gaussian 
white noise. In principle, the noise components could be correlated. But since the stretch 
current is itself conserved, its fluctuations will be taken care of by the momentum equation. 
Momentum and energy currents have different signature under time reversal, hence their 
cross correlation vanishes. As a result, there is a fluctuating momentum current of strength 
(Ju and an independent energy current of strength a^- According to Onsager, noise is linked 
to dissipation as modeled by a diffusive term. Thus the linearized equations (2.27) are 
extended to 

dtu{x,t) + dx{Au{x,t) — dxDu{x,t) + B^{x,t)) =0. (3.1) 

Here ^{x,t) is standard white noise with covariance 

= 6aa’S{x - x')6{t -t') (3.2) 

and, as argued, the noise strength matrix is diagonal as 

H = diag(0, cTu, cTe). (3.3) 

To distinguish the linearized Euler equations (2.27) from the Langevin equations (3.1), 
we use u = (ui, U 2 , U 3 ) for the fluctuating helds. 

From the introduction, we know already that a Gaussian fluctuation theory will fail. 
But still, it is useful to hrst explore the structure of the Langevin equation (3.1). The 
stationary measures for (3.1) are spatial white noise with arbitrary mean. Since small 
deviations from uniformity are considered, we always impose mean zero. Then the com¬ 
ponents are correlated as 


(u„(x)Ua/(x')) = Caa'Kx - x') . (3.4) 

Stationarity relates the linear drift and the noise strength through the steady state co- 
variance as 

- (AC - CA^)dx + (DC + CD^)d^ = BB^dl. (3.5) 

The hrst term vanishes by (2.30) and the diffusion matrix is uniquely determined as 

/ 0 0 0 \ 

D = 0 D, 0 . (3.6) 

V^e 0 dJ 

with = —(ro; Vb)p,/ 3 (ro; ro)p^Ile- Here > 0 is the momentum and He > 0 the 
energy diffusion coefficient, which are related to the noise strength as 

'^u = {po',Po)p,/3Du , cTg = (eo; eo)p,/3He. (3.7) 

We still have to establish that the stationary measure (3.4) is unique and is approached 
in the limit t —)• oo. For this purpose it suffices that the 3x3 matrix YlirkA — (27rfc)^H has 
its eigenvalues in the left hand complex plane, where for convenience we have switched 
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to Fourier space with respect to x. If one drops D, then \2TTkA has the eigenval¬ 
ues i27rfcc(—1, 0,1). Hence one can use hrst order perturbation theory with respect to 
— {27ikYD, which is given by D'll’a), where = Cai^a and AAt/ja = Cq'^q, are the 
right and left eigenvectors of A, as listed in [29]. One simply has to follow the dehnitions 
and express every term through the cumulants of tq and Vq. As to be expected, the ma¬ 
trix elements from above have a dehnite value and the eigenvalues are shifted into the left 
hand complex plane. Similarly, —D has eigenvalues 0, — —D^ and the zero eigenvalue 
is shifted to the left by second order perturbation in i(27r/c)“^A. The only condition is 
the strict positivity of Du,De. 

Based on (3.1) one computes the stationary space-time covariance, which most easily 
is written in Fourier space, 

Saa'{x,t) = Mx,tWi0,0)) = j . (3.8) 

To extract the long time behavior it is convenient to transform to normal modes. But 
before, we have to introduce a more systematic notation. We will use the superscript ** 
for a normal mode quantity. Thus for the anharmonic chain 

S\j,t) = RS{j,t)R^, = {{RgUj,t){Rg)aiO,0))p,f,. (3.9) 

The hydrodynamic fluctuation helds are dehned on the continuum, thus functions of x, t, 
and we write 


*S'aa'(a;, f) = (ua(a;, f)uQ,/(0, 0)) , S^{x,t) = RS{x,t)R^ ■ (3.10) 


Correspondingly A^ = RAR~^ = diag(—c, 0, c), D‘^ = RDR~^, B'^ = RB. Note that 
u{x, t) will change its meaning when switching from linear to nonlinear fluctuating hydro¬ 
dynamics. 

In normal mode representation Eq. (3.8) becomes 



^i27Tkx^—it27TkA^ — \t\{27Tk)‘^ 


(3.11) 


The leading term, it27rkA\ is diagonal, while the diffusion matrix D'^ couples the compo¬ 
nents. But for large t the peaks are far apart and the cross terms become small. More 
formally we split D’^ = D^ia + Dos and regard the off-diagonal part Dos as perturbation. 
When expanding, one notes that the off-diagonal terms carry an oscillating factor with 
frequency — Cq,/, a ^ a'. Hence these terms decay quickly and 




i2nkx^—it2TTkcc, — \t\ {2nk)'^ 


(3.12) 


for large t. Each peak has a Gaussian shape function which broadens as ■ 

Besides the peak structure, we have gained a second important insight. Since the 
peaks travel with distinct velocities, on the linearized level the three-component system 
decouples into three scalar equations, provided it is written in normal modes. 
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The linear fluctuation theory should be tested against adding nonlinear terms. The 
computation of the Introduction indicates that in one dimension the quadratic part of 
the Euler current is a relevant perturbation. This can be seen even more directly by 
rescaling u(x,t) to large space-time scales as Ue(x, t) = e“^Ue(e“^a;, and counting the 
powers of the nonlinear terms. To have the correct t = 0 covariance requires b = 1/2. 
The quadratic terms of the Euler currents are relevant, yielding 2 ; = 3/2, and the scaling 
exponents of cubic terms are only marginally relevant, while a possible dependence of 
D, B on u{x,t) can be ignored. Thus, we retain the linear part (3.1) but expand the Euler 
currents including second order in u, which turns (3.1) into the equations of nonlinear 
fluctuating hydrodynamics, 

^tUi - d^U2 = 0 , 



To explore their consequences is a more demanding task than solving the linear Langevin 
equation and the results of the analysis will be more fragmentary. 

Stationary measure for the physical fields. Adding quadratic terms could change 
drastically the character of the solution. To And out we first attempt to investigate the 
time-stationary measure. The vector held of the nonlinear part of (3.13) reads 



Formally the drift is divergence free, since 

J = deP j dx{-d^U2 + dxU2) = 0 (3.15) 

and the infinite dimensional Lebesgue measure is invariant under the flow generated by 
F. Since the equilibrium measure is a product, a natural ansatz for the invariant measure 
is Gaussian white noise retaining the physical susceptibility. 



(3.16) 


As established before, for the linear Langevin equation this measure is stationary. Thus, 
to And out whether it is also stationary for (3.13), we only have to check the invariance 
under the nonlinear flow 


d 

dt 



3 


3 


/ 


dxaoU2dx{id^P)ul - {deP)ul + {d‘^P)ul -F 2{dideP)uiU3) 



(3.17) 
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where Oq = (C'“^) 22 , cl 2 = <23 = (C'~^) 33 - The term cubic in U 2 vanishes by the 

same argument as for the one-component case. All other terms are linear in U 2 , thus their 
sum has to vanish point-wise, 


0 = {aodjP - a2diP)d:^ul + {aod'^P - a3deP)d:^ul 

+2{ao{dideP)d^{uiU3) - a2{deP)u3d^ui - a3{diP)uid^U3) . (3.18) 


This leads to the constraints on the coefficients as 

aodJP = a2diP , aQd‘^P = Os^e-P , aodid^P = , aodid^P = a 2 pe-P • (3.19) 

There are four constraints and hve partial derivatives which may be regarded as inde¬ 
pendent parameters. Thus one would expect that there is a sub-manifold in V, P, space 
for which (3.19) can be satished. We will come back to another representation of these 
constraints below. Away from the special subset of invariant Gaussian measures, we have 
no tools, but one would hope that the invariant measure still has a hnite correlation length 
and exponential mixing. Based on the mechanical model, it is suggestive to assume U 2 
to be independent of ui, U 3 and to have white noise statistics. But this forces again the 
constraints (3.19) and results in the same Gaussian measure as before. 

A basic property of the mechanical model is invariance under time reversal. On the 
level of fluctuating hydrodynamics this translates to the following property: We £x a time 
window [0,T]. Then, in the stationary process, the trajectories 


(ui(t), U 2 (t), U 3 (t)) and (ui(T - t), -U 2 (T - t), U 3 (T - t)) 


(3.20) 


with 0 < t < T have the same probability. To check (3.20) requires some information 
on the invariant measure. But under the Gaussian measure (3.16), hence assuming the 
validity of the constraints, it can be shown that time reversal invariance indeed holds. 

Transformation to normal modes. To proceed further, it is convenient to write (3.13) 
in vector form. 


dtu{x, t) -I- dx{Au{x, t) + ^(u, Hu) — dxDu{x, t) + B^{x, t)) = 0 , (3.21) 


where H is the vector consisting of the Hessians of the currents with derivatives evaluated 
at the background values (£, 0 , e). 


3 




(3.22) 


As for the linear Langevin equation we transform to normal modes through 


cj) = Ru . 


(3.23) 


Then 


dt<Pa + dx[Ca4>a + (0, G“0) — dx{D'^(j))a “t" (P^^)a) — 0 . 


(3.24) 
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By construction = 2DK The nonlinear coupling constants, denoted by G, are 

defined by 

3 

G" = i ^ (3.25) 

with the notation R~"^ = 

Since derived from a chain, the couplings are not completely arbitrary, but satisfy the 
symmetries 


GL 

— ^7/3 : : ^-10 ~ ^01 5 

= —G° , G°^ = 0 otherwise . 

(3.26) 

In particular note that 

o o 
o 

o 

(3.27) 


always, while G\i = are generically different from 0. This property signals that 

the heat peak will behave differently from the sound peaks. The G-couplings are listed 
in [29] and as a function of P, /3 expressed in cumulants up to third order in ro, Vq. The 
algebra is somewhat messy. But there is a short MATHEMATICA program available [32] 
which, for given P, /?, V, computes all coupling constants, including the matrices G, A, R. 

We return to the issue of Gaussian time-stationary measures, where we regard the 
coefficients G, P**, B^ as arbitrary, up to 2P** = B^B'A^, ignoring for a while their particular 
origin. The Langevin equation (3.24) is slightly formal. To have a well-defined evolution, 
we discretize space by a lattice of N sites. The field (j){x,t) then becomes with 

components 0j^a(f), j = a = 1,2,3. The spatial finite difference operator is 

denoted by dj, djfj = fj+i — fj, with transpose djfj = fj-i — fj. Then the discretized 
equations of fiuctuating hydrodynamics read 

dt(l>j,a + dj{Ca4’j,a + ~ ^ (3.28) 

with (pj = (j)N+ji ^0 = ^N, where are independent Gaussian white noises with covari¬ 
ance 

= Sjfdaa'S{t - f'). (3.29) 

The diffusion matrix and noise strength B^ act on components, while the difference 
operator dj acts on the lattice site index j. 

is quadratic in 0. But let us first consider the case A/)> = 0. Then ((>j,a{t) is a 
Gaussian process. The noise strength has been chosen such that one invariant measure is 
the Gaussian 

N 3 N 3 

n n =pg( 0) n n • (^-so) 

j=l a=l j=l a=l 

Because of the conservation laws, the hyperplanes 

N 

= Npc (3.31) 

i=i 
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are invariant and on each hyperplane there is a Gaussian process with a unique invariant 
measure given by (3.30) conditioned on that hyperplane. For large N it would become 
independent Gaussians with mean our interest being the case of zero mean, = 0. 
The generator of the diffusion process (3.28) with Mj^a = 0 is given by 

N 3 3 

j=l 0=1 a,Q'=l 

The invariance of Pg{(P) can be checked through 

L*pg(0) = O, (3.33) 

where * is the adjoint with respect to the flat volume measure. Furthermore linear func¬ 
tions evolve to linear functions according to 

N 3 

e^°Vi,a = , (3.34) 

j^=l 0-'=! 


where the matrix A = —dj ® diag(ci, C 2 , C 3 ) — djdj ® the hrst factor acting on j and 
the second on a. 

We now add the nonlinearity In general, this will modify the time-stationary 

measure and we have little control how. Therefore we propose to choose A/)_q, such that 
the corresponding vector held djNj^a is divergence free [9]. If Afj^a depends only on the 
held at sites j and j -l- 1 , then the unique solution reads 


3 

A/j^Q. = 2 'y ^ -1- T 0j+l^2'0i+l,7') • 

7 , 7'=1 


(3.35) 


For pg to be left invariant under the deterministic how generated by the vector held —djM 
requires 

N 3 

FiPg = 0 , Li = — (3.36) 

i=i «=i 


which implies 


N 3 

i=i «=i 


(3.37) 


and thus the constraints 

GJ, = Gip = G”^) (3.38) 

for all a,f3,'j = 1,2, 3, where in brackets we added the symmetry which holds by dehnition. 
Denoting the generator of the Langevin equation (3.28) by 


L — Lq + Li, 


(3.39) 
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one concludes L*pQ = 0, i.e. the time-invariance of pq. 
In the continuum limit the condition (3.37) reads 


^ f dx(j)a{x)d^{(j)0{x)(j)^{x)) = 0 , 

«./3i7=1 

where By partial integration 

3 „ 3 


(3.40) 


2 ^ J dx(paix)(j)i3ix)da:(l)^{x) = - ^ G^^ J dx(j)fs{x)(j)^{x)d^(paix) ■ (3.41) 


a,/3,7=l 


a,/3,7=l 


Hence (3.40) is satished only if = G'^^, which is the condition (3.38) claimed for the 
discrete setting. 

(3.38) is the generalization of (3.19), which is specihc for the anharmonic chain. In 
fact, while abstractly true, it is not so easy to verify directly. But now we can argue more 
convincingly why one should be allowed to continue with assuming the validity of the 
constraints (3.38). As we will discuss in the next section, the leading coupling constants 
are of the form while the sub-leading couplings have equal lower indices, G^^, 7 7^ «• 
The off-diagonal matrix elements are irrelevant for the large scale behavior. When one 
does the counting, all leading and sub-leading couplings can be chosen freely and the 
irrelevant couplings can be adjusted so that the constraint (3.38) is satished. Appealing 
to universality, the large space-time behavior should not depend on that particular choice. 
We expect that for general G the true time-stationary measure will have short range 
correlations and nonlinear huctuating hydrodynamics remains a valid approximation to 
the dynamics of the anharmonic chain. 

In related problem settings, a different point of view has been suggested [33, 34]. 
Firstly one notes that the Gaussian stationary measure (3.16), hence also (3.30), is simply 
inherited from the canonical equilibrium measure. In this respect there is no choice. 
Also the nonlinear Euler currents are on the safe side. But the remaining terms are 
phenomenological to some extent. D^, B'^ could depend on 0 itself. One could also include 
higher derivative terms. In fact, one could try to choose the nonlinearities precisely in 
such a way that the dynamics is invariant under time-reversal and leaves the Gaussian 
measure invariant. The program as such may be easily endorsed, but so far I have not 
seen a convincing handling of the details. 


4 Mode-coupling theory 

Decoupling hypothesis. For the linear equations the normal modes decouple for long 
times. The hypothesis claims that such property persists when adding the quadratic 
nonlinearities. For the precise phrasing, we have to be somewhat careful. We consider a 
hxed component, a, in normal mode representation. It travels with velocity Cq,, which is 
assumed to be distinct from all other mode velocities. If G°^ ^ 0, then for the purpose of 
computing correlations of mode a at large scales, one can use the scalar conservation law 

dtK + = 0 , (4.1) 
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which coincides with the stochastic Burgers equation (1.19). If decoupling holds, one 
has the exact asymptotics as stated in ( 1 . 21 ) with A = 2 a/ 2|G"„|. The universal scaling 
function /kpz is tabulated in [35], denoted there by /. /kpz > 0, / dxfKPz{x) = 1, 
fKPz{x) = fKPz{—x), f dxfKPz(x)x^ ~ 0.510523. /kpz looks like a Gaussian with a large 
|a;| decay as exp[—0.295|a;|^]. Plots are provided in [17, 35]. 

For an anharmonic chain, Gqq = 0 always and the decoupling hypothesis applies only 
to the sound peaks, provided G\^ = 7 ^ 0 which generically is the case. If 7 ^ 0, 

then the exact scaling form is 

St{x,t) = {Kt)-^/^fKPz{{Kt)-^/\x - act)) , As = 2 ^|G^J , (4.2) 

a = ±1. To hnd out about the scaling behavior of the heat mode other methods have to 
be developed. 

For one-dimensional fluids, van Beijeren [36] follows the scheme developed in [ 6 ] and 
arrived hrst at the prediction (4.2) together with the Levy 5/3 heat peak to be discussed 
below. In [36] no Langevin equations appear. I regard them as a useful intermediate 
step valid on a mesoscopic scale. In the Langevin form the theory can be applied to a 
large class of one-dimensional systems. As a tool, fluctuating hydrodynamics has been 
proposed considerably earlier [37] and used to predict the decay of the total energy 
current correlation. 

One-loop, diagonal, and small overlap approximations. We return to the Langevin 
equation (3.28) and consider the mean zero, stationary </yQ-(t) process with pq as t = 0 
measure. The stationary covariance reads 

>S'L'(T^) = (‘/’j>(^)0O,a'(O)) = (0o,a'e^Vya)eq, t>0. (4.3) 

On the left, (■) denotes the average with respect to the stationary 0j,a(^) process and on 
the right (■)eq refers to the average with respect to pc- By construction 

~ ■ ( 4 - 4 ) 

The time derivative reads 

= (0O,a'(e“Lo0ya))eq + (0O,a'(e^*Li0j-„))eq . (4.5) 

We insert ^ 

f (4.6) 

Jo 

in the second summand of (4.5). The term containing only e^°* is cubic in the time zero 
helds and hence its average vanishes. Therefore one arrives at 

(e^^l/j-«))eq . (4.7) 

For the adjoint of ^gg ( 3 . 34 ) and for the adjoint of Li we use 

„LlF(0))eq = -((Li0j-o)F(0))eq , (4.8) 
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(4.9) 


which both rely on (•)eq being the average with respect to pq. Furthermore 
Inserting in (4.7) one arrives at the identity 

• (4-10) 

To obtain a closed equation for S'** we note that the average 

{dj>J\ff^a'e^‘’djJ\fj^a)eq = {djAfj^a{s)dfJ\ff ^a'{0)) (4.11) 

is a four-point correlation. We invoke the Gaussian factorization as 

(0(s)(/>(s)(/>(O)0(O)) = ((^(s)0(s))(0(O)0(O)) + 2(0(s)(/)(O))(0(s)0(O)). (4.12) 


The hrst summand vanishes because of the difference operator dj. Secondly we replace 
the bare propagator by the interacting propagator S**(f — s), which corresponds 

to a partial resummation of the perturbation series in G. Finally we take a limit of 
zero lattice spacing. This step could be avoided, and is done so in our numerical scheme 
for the mode-coupling equations. We could also maintain the ring geometry which, for 
example, would allow to investigate collisions between the moving peaks. Universality is 
only expected for large j, t, hence in the limit of zero lattice spacing. The continuum limit 
of S'**(j,f) is denoted by S^{x,t), x G M. With these steps we arrive at the mode-coupling 
equation 


3 

(( - cJaa'dx + D^^,dl)sl^{x,t) 

a'=l 

+ [ ds [ dyM^a'{y,s)dlSl,^{x - y,t - s)) (4.13) 

Jo Jr ' 

with the memory kernel 

3 

M„„/(x,t)=2 ^ G^,yG^,',y,S'J,^„(x,t)S'“,y,(x,f). (4.14) 


In numerical simulations of both, the mechanical model of anharmonic chains and the 
mode-coupling equations, it is consistently observed that (j, t) becomes approximately 
diagonal fairly rapidly. To analyse the long time asymptotics on the basis of (4.13) we 
therefore rely on the diagonal approximation 

~ Saa’fa{x,t) . (4.15) 


Then fa{x,0) = 6{x) and the f^s satisfy 

dtfa{x,t) = {-Cadx + Dl^dl)fa{x,t) + [ ds [ dydlfa{x-y,t-s)Maa{y,s), (4.16) 
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a = — 1 , 0 , 1 , with memory kernel 


M^^{x,t) = 2 {G^yYf^{x,t)fY{x,t). (4.17) 

7,7'=0,±1 

The solution to (4.16) has two sound peaks centered at ±ct and the heat peak sitting 
at 0. All three peaks have a width much less than ct. But then, in case 7 7 ^ 7 ', the 
product f^{x,t)fY{x,t) ~ 0 for large t. Hence for the memory kernel (4.17) we invoke a 
small overlap approximation as 

Maa{x, t) ~ t) = 2 Y tf , (4.18) 

7 = 0,±1 


which is to be inserted in Eq. (4.16). 

Numerical simulations of the mode-coupling equations. When starting this project 
together with Christian Mendl, in the summer of 2012 we spent many days in numerically 
simulating the mode coupling equations with initial conditions *S'^c^/(j,0) = Saa'^oj- Only 
a few plots are in print [38], simply because there is such a large parameter space and it 
is not clear where to start and where to end. Still, for our own understanding this period 
was extremely helpful. Mostly we simulated in Fourier space. System size was up to 400. 
Speeds were of order 1 , thereby limiting the simulation time to about 200, the time of the 
hrst peak collision. For such sizes the simulations are fast and many variations could be 
explored. We started from the scalar equation, to be discussed below, moved up to two 
modes, and eventually to three modes with parameters taken from an actual a n h armonic 
chain. | was either 0 or somewhere in the range 0.3 to 2.5. D'^ is a free parameter 

which was varied from 0 to |G",„,|/2. We always simulated the complete matrix-valued 
mode-coupling equations (4.13). Our main Endings can be summarized as: 

(i) For a large range of parameters, the diagonal approximation in generally failed for 
short times, but was quickly restored with the oE-diagonal elements being at most 10% 
of the diagonal ones. 

The results were fairly insensitive to the choice of DK In fact, = 0 works also well. 
Apparently the memory term generates already enough dissipation. 

(Hi) We varied the overlap coefficients with 7 % 7 '. Over the time scale of the 
simulation no substantial changes were observed. 

All these findings confirm the approximations proposed. 

As our biggest surprise, except for trivial cases we never reached the asymptotic regime. 
The peak structure develops fairly rapidly. The peak shape then changes slowly, roughly 
consistent with the predicted scaling exponents, but it does not reach a self-similar form. 
For example in the case of the sound peak, on the scale rather than being symmetric, 
as claimed by (4.2), it is still badly distorted, tilted away from the central peak with rapid 
decay outside the sound cone but rather slow power law type of decay towards the central 
peak. To improve one would have to simulate larger system sizes and longer times. But 
then numerical simulations become heavy and the fun evaporates. More attention can be 
achieved by molecular dynamics (MD) simulation of the mechanical chain. 
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For given parameters V, P, (3 one easily computes all the required coefficients. So one 
goal was to run a MD and put the results on top of the ones from a simulation of the mode 
coupling equations. For this to be a reasonable program, one would have to simulate the 
mode-coupling equations for sizes of iV = 4000 and more, which we never attempted. 

For the scalar case the situation is much simpler. The mode-coupling equation takes 
the form 

dtf{x,t) = Ddlf{x,t)+ 20“^ [ ds [ dydlf{x-y,t-s)f{y,sf , (4.19) 

Jo Jr 

which is the one-loop approximation for the stochastic Burgers equation [39] . For large 
x,t, its solution with initial condition f{x,0) = 6{x) takes the scaling form 

f{x,t) = . (4.20) 

Inserting in (4.19), one hrst hnds the non-universal scaling coefficient 

K = 2V2\G:J. ( 4 . 21 ) 

Secondly /me, the Fourier transform of /me, is dehned as solution of the hxed point equation 

IfLM = [ ds/mc((l - s)^/^^) [ dqfmc{s^^^{w - g))/mc(s^^^g) (4.22) 

Jo Jr 

with w > 0 and /mc(0) = 1, /me(0) = 0- 

Eq. (4.22) is based on the closure assumption (4.12) and there is no reason to infer 
that it is exact. However from our numerical simulations we conclude that /kpz differs 
from /me by a few percent only. We regard the scalar case as a strong support for the 
entire approach. But for several components the large hnite size effects prohibit one to 
arrive at a similarly simple claim. 

Asymptotic self-similarity. Within mode-coupling the asymptotic shape function for 
the sound peaks is given by 

fa{x,t) = fmc{iKt)~‘^^^{x - (TCt)) , (4.23) 

a = ±1. For the heat peak we employ (4.16) together with (4.18), using as an input that 
the asymptotic form of is known already. In fact the scaling exponent for /o- is crucial, 
but the precise shape of enters only mildly. Hence, again switching to Fourier space, 
one has to solve 

dtfo{k,t) = -Dl{2TTkffoik,t) 

f dsfo{k,t-s) f dqfa{k-q,s)fa{q,s), (4.24) 
a=±i do Jr 

/o(fc,0) = 1. For /o- one inserts the asymptotic result (4.20). (4.24) is a linear equation 
which is solved through Laplace transform with the result 

/„(fc,i) (4.25) 


26 


where 


Ah = Ag f dtt cos(27rct) f dxf^^ixf 

Jo Jm. 

= As“^''^(G'L)^(4vr)^(27rc)-^/^i7r-^—^ f dxf^^{xf (4.26) 

i COSl^gJ Jr 

and we used the symmetry see [29] for details. (4.25) is the Fourier 

transform of the symmetric a-stable distribution with exponent a = 5/3, also known as 
Levy distribution. In real space the asymptotics reads, for |x| > (Aht)^^®, 

/o(a;, t) ~ 7r“^Aht|a;|“®/^ . (4.27) 

/me is a smooth function with rapid decay. On the other hand, /o has fat tails and 
its variance is divergent. According to (4.27), at a; = act the heat peak /o((Tcf,f) = 
This explains why there is still coupling between /q and /o-, despite the 
large spatial separation. In fact, numerically one observes that beyond the sound cone, 
X = ±ct, the solution decays exponentially fast. As t becomes large the tails of /o are 
build up between the two sound peaks, so to speak they unveil the Levy distribution. 

In (4.24) we could also insert for f^. the exact scaling function /kpz, which would 
slightly modify Ah. fhevy is an approximation, just as /me- But the MD simulations 
display so convincingly the Levy distribution that one might be willing to regard it as 
exact. If so, the exact Ah must be based on /kpz- To obtain the correlations of the 
physical helds, one has to use 


Sij,t) = R-^S\j,t)R-^. (4.28) 

In particular the correlations of the physical helds are given through 

SUj,t)= (4.29) 

( 7 = 0,±1 

where for /o- the asymptotic scaling form is inserted. Then asymptotically the i-i and 
the e-e correlations show generically all three peaks. However, for the u-u correlations the 
central peak is missing asymptotically, since (i?~^) 2 o = 0 independently of the interaction 
potential V. 

We note that the coefficient D'^ does not appear in the asymptotic scaling form, of 
course neither = 2DK This result is consistent with the picture that noise and dis¬ 

sipation are required to maintain the correct local stationary measure with susceptibility 
C. The long time asymptotics is however governed by the nonlinearities. 

No signal beyond the sound cone. Physically the sound speed is an upper limit 
for the propagation of small disturbances. Since the initial state has a hnite correlation 
length, one would expect that towards the outside of the sound cone correlations decay 
exponentially, while inside the sound cone there seems to be no particular restriction. As 
a consistency check, one would hope that such a general feature is properly reproduced 
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by mode-coupling. Their numerical solutions conform with this expectation, at least for 
the small system sizes explored. But for the scaling limit one has to let t —)■ oo and the 
decay information seems to be lost. However there is still a somewhat subtle trace. 

To explain, we have to hrst recall some properties of Levy stable distributions. Except 
for trivial rescalings, they are characterized by two parameters, traditionally called a, (3, 
where for simplicity we momentarily stick to this convention, without too much risk of 
confusion. The probability density has a simple form in Fourier space, 

/Levy,Q,/ 3 (fc) = exp ( - - i/3 tan(| 7 ra)sgn(/c)]). (4.30) 

The parameter a controls the steepness, 0 < a < 2, while (3 controls the asymmetry, 
1/31 < 1. For |/3| > 1 the Fourier integral no longer defines a non-negative function. At the 
singular point a = 2 only /3 = 0 is admitted and the probability density is a Gaussian. If 
|/3| < 1, the asymptotic decay of fhevy,a,i 3 {x) is determined by a and is given by 
for |a;| oo. At |/3| = 1 the two tails show different decay. The functions corresponding 
to /3 = 1 and /3 = —1 are mirror images, for /3 = 1 the slow decay being for x —)■ —oo 
and still as For 0 < a < 1, fLevy,a,i{x) = 0 for x > 0, while for 1 < a < 2 the 

decay becomes stretched exponential as exp(—with known constant Cq. We 
refer to [40] for more details. 

For the heat peak we obtained the symmetric Levy distribution because the sound 
peaks are reflection symmetric, implying ci = —c_i and If hypo¬ 

thetically we would choose distinct couplings, or ci 7 ^ —c_i, then this imbalance would 
produce a /3 7 ^ 0. If one of the sound peaks would be completely missing, as in the case 
for a system with only two conserved helds, then necessarily |/3| = 1. In accordance with 
the physical principle, the sign of /3 is such that the fast decay of /Levy,a,±i(a^) is towards 
the outside of the sound cone, while the slow decay is towards the single sound peak. For 
hnite t, this slow decay will be cut by the sound peak. Thus the scaling solution of the 
mode-coupling equations reproduces the rapid decay towards the exterior of the sound 
cone. This is a completely general fact, any number of components and any G [41]. 

Dynamical phase diagram. As already indicated through the particular case Gqq = 0, 
the large scale structure of the solution depends on whether = 0 or not. One extreme 
case would be G"„ 7 ^ 0 for all a, implying that the three peaks have KPZ scaling be¬ 
havior. The other extreme is = 0 for all a, 7 , resulting in all peaks to have diffusive 
broadening. For the case of only two modes, the full phase diagram has seven distinct 
phases, with unexpected details worked out in [42]. For the general case of n components 
the long time asymptotics is completely classihed in [43, 41]. Anharmonic chains have 
special symmetries and not all possible couplings G can be realized. Given that Gqq = 0 
and because the sound peaks are symmetric, to have a distinct scaling requires 

G}i=0, (4.31) 

which can be realized. The behavior is then determined by the value of the remaining 
diagonal matrix elements. For the central peak one finds 

^GL > 0, 
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(4.32) 


while for the sound peak diagonals, Gii_i = —Gi/, Gqo = there seems to be no 

particular restriction. In principle, there could be sort of accidental zeros of which 
are then difficult to locate. A more direct approach starts from the observation that the G 
coefficients are expressed through cumulants in tq, Vq . If the integrands are antisymmetric 
under reflection, many terms vanish. The precise condition on the potential is to have 
some oo, Pq such that 


V{x — ao) + Pqx = V {—X — ao) — Pqx (4.33) 

for all X. Then for arbitrary (3 and P = Po, one hnds 

= 0 , = 0 , Gqq = —Gqq^ = 0, (4.34) 

while > 0. The standard examples for (4.33) to hold are the FPU chain with no 

cubic interaction term, the /9-chain, and the square well potential with alternating masses, 
both at zero pressure. 

Under (4.34) the heat mode is coupled to the sound mode, but there is no back reaction 
from the sound mode. Hence the sound peak is diffusive with scaling function 

UX, t) = ^ 

Ds is a transport coefficient. It can dehned through a Green-Kubo formula, which also 
means that no reasonably explicit answer can be expected. The feed back of the sound 
peak to the central peak follows by the same computation as before, with the result 

/o(fc,t)=e-l'=l'''"‘'S (4.36) 

where 

poo 

= cos{t){2^/^T)~^. (4.37) 

Jo 

Since 3/2 < 5/3, the density fo{x,t) turns out to be broader than the Levy 5/3 from the 
dynamical phase with G\i ^ 0. 

In testing nonlinear fluctuating hydrodynamics almost subconsciously one tries to 
conhrm (or not) the scaling exponents, resp. functions. This can be difficult because 
of limited size. The dynamical phase diagram offers a different option. For exceptional 
points in the phase diagram, without too precise a verihcation of the scaling, one should 
hnd that the standard scaling exponent does not properly £t the data. Such qualitative 
property is possibly more easy to access. In Fig. 2 we display heat and sound peak for a 
FPU chain with Gj^ = 0. 

5 Molecular dynamics simulations 

In 1953 Fermi, Pasta and Ulam, technically supported by Tsingou, simulated 32 particles 
interacting through a quartic potential at the extremely low energy of e = 5 x 10“^ per 
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Figure 2: Scaling plot of heat and sound peak for a FPU chain with N = 8192, potential 
parameters a = 0, (3 = 1, pressure P = I, and temperature = 1 . 

particle (above the ground state energy) [25]. They injected energy in the highest Fourier 
mode and were looking for equipartition of the modes at long times. However they found 
quasi-periodic motion with time averages settling to some dehnite value different from 
equipartition. The observed quasi-periodicity triggered the connection to KAM tori, the 
discovery of integrable systems with many degrees of freedom, and to the development 
of the theory of solitons and breathers [44, 45]. For sure, a rich harvest, see [46] for a 
historical perspective. Later on Izrailev and Chirikov [47] repeated the simulation at the 
higher energy e = 5 x 10“^ and observed equipartition. 

Anomalous transport surfaced much later [48], see the reviews [49, 50]. One connects 
the two ends of the chain to thermal reservoirs. To explain, the end particles are tied 
down as go = 0 , gw+i = 0 , and to the equations of motion for the boundary particles one 
adds Langevin terms as 

<?i = -7Pi + , gjv = -7Pv + ^27^+^+ > (5-1) 

where 7 is a friction constant, T± are the boundary temperatures, and ^±{t) are indepen¬ 
dent standard Gaussian white noises. For T_ = T+ the system settles in the canonical 
equilibrium state. But for T_ 7 ^ T+ there is a non-trivial steady state with a non-zero 
energy flux je{N) depending on the length, N, of the chain. For regular heat transport, 
Fourier’s law implies je{N) — CoN~^. However for FPU chains one hnds an enhanced 
transport as 

je(A^)^coiV-i+“(T_-T+) (5.2) 

with an exponent a characterizing the anomaly. [A further a, but better to stick to 
standard conventions]. Over the last two decades many MD simulations have been imple¬ 
mented for a wide variety of one-dimensional systems. Early results indicated a = 2/5, 
but since about 2003 the common evidence pointed towards a = 1/3 or at least close to 

it. 


30 

















Nonlinear fluctuating dynamics can also deal with such open chains, at least in prin¬ 
ciple. One would impose energy imbalance boundary conditions as -^(0,^) = 0 = £{L,t), 
u(0,t) = 0 = u(L,t), but c(0,t) = c_ and t{L,t) = c+ and tries to investigate the steady 
state. Unfortunately, at least for the moment, we have no powerful techniques to deal 
with this problem. On the other hand it is argued [49] that the energy flux is related to 
a Green-Kubo formula by 

pN/c 

UN) ~ / dt{{JU); Jsm - (J3(oo); J3(0))) . (5.3) 

Jo 

Under the time integral appears the total energy current correlation in thermal equilibrium 
with its possibly non-zero value at t = oo subtracted, see Section 6 for more explanations. 
The decay of such correlation can be predicted by mode-coupling. In fact, we will conhrm 
the value a = 1/3. But the argument is subtle because it is only indirectly related to the 
spreading of the heat peak which is on scale 

Because of the relation (5.3), in many MD simulations the total energy current corre¬ 
lation ( 573 (f);^ 73 ( 0 )) is measured as an addition to steady state transport [51, 52]. There 
are also MD simulations exclusively focussed on momentum and energy current correla¬ 
tions [53]. However simulations of correlations of the conserved helds have been fairly 
scarce until recently. The peak structure was noted already in [54], see also [55, 56]. But 
surprisingly enough, even such a basic issue as the quantitative comparison between the 
measured speed of sound and formula (2.31) is apparently not a routine check. So far 
there have been three independent sets of simulations with the specific aim to check the 
predictions from mode-coupling. For the details the reader is encouraged to look at the 
original papers. I will try to compare so to reach some sort of conclusion. 

The first and second set are FPU chains with either symmetric or asymmetric poten¬ 
tial, both the a[3 and the pure fJ chain [57, 58]. In this case one has to integrate numer¬ 
ically the differential equations governing the evolution, for which both a velocity-Verlet 
algorithm and a fourth order symplectic Runge-Kutta algorithm are used. The third set 
consists of chains with a piecewise constant potential [59]. We call them hard-collision, 
since the force is zero except for 5-spikes and the dynamics proceeds from collision to colli¬ 
sion. Now one has to develop an efficient algorithm by which one hnds the time-wise next 
collision. Except for rounding, there is no discretization of time. The simplest example 
is the hard-point potential, VU^) = C )0 for x < 0 and VU^) = 0 for x > 0. A variant 
is the infinite square well potential, Uw(a^) = C )0 for x < 0, x > a and VU^) = 0 
0 < X < a [60, 59]. In this case two neighboring particles at the maximal distance a are 
reflected inwards as if connected by a massless string of hnite length a. For both models 
the dynamics remembers the initial velocities. To have only the standard conservation 
laws, one imposes alternating masses, say mj = mo for even j and mj = mi for odd j. In 
both models the unit cell then contains two particles and the scheme explained before has 
to be extended. But at the very end the difference is minimal. A further variant is the 
shoulder potential Uh(3^), for which Uh(3^) = cxd for |x| < |, Uh(3^) = for | < |a^| < 1, 
and I4h(x) = 0 for |x| > 1. The potential is either repulsive, eo > 0, or attractive eo < 0. 
Exploratory studies of the latter case indicate that the convergence is slower than for 
the extensivlely studied attractive case. Particles interacting with such a potential can 
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be viewed also as a hard-core fluid with a short range potential part and thus serves as 
a bridge between one-dimensional fluids and anharmonic solids. The collisions resulting 
from the potential step make the model non-integrable. 

FPU chains with an even potential at P = 0 constitute a distinct dynamical phase. 
Such phase is absent for the hard-point and the square shoulder potentials. But the square 
well at P = 0 has the same properties as can be seen from taking oq = a/2 in (4.33). 

Current system sizes are N = 2^^ to 2^^, even size 2^® = 65, 536 has been attempted 
[61]. Periodic boundary conditions are imposed. Time is restricted to t < tmax = N/2c, 
the time of the first collision between sound peaks. For given potential, one has to decide 
on the thermodynamic parameters, P, (3. One constraint is to have c approximately in the 
range 1... 2 in order to have a sufficiently long time span available. Secondly one would 
like to be well away from integrability. This leads to an energy per particle of order 1 in 
the models from above. A related issue are the coupling matrices G, which should not 
be too small, at least for the relevant couplings. For the models from above they are 
tabulated. The relevant couplings show quite some variability taking values in the range 
0.1... 3.4. At the very end, one has to make a physically reasonable choice, perhaps use 
the same parameters as previous MD simulations so to have the possibility to compare. A 
systematic study of the dependence on P, (3 seems to be too costly and most likely not so 
interesting. But it does make sense to probe values at the border. For some of the FPU 
simulations the energy per particle is chosen as e = 0.1, at which point nonlinearities are 
small [57]. There are also very extensive simulations at the even lower energy e = 5 x 10““^, 
with the goal to explore the route to equipartition, which is a somewhat distinct issue 
162 ], 

Once all parameters are fixed, there are several options to run the simulation. Since 
in canonical equilibrium {rj,pj,j = 1,..., A^} are independent random variables, one can 
sample the initial conditions through a random number generator. For the hard-collision 
potentials the geometric constraints are still simple enough for allowing one to generate 
the microcanonical ensemble by Monte Carlo methods. In our simulations the correlator 
S hardly depends on the choice of the ensemble. With such generated random initial 
data the equations of motion are simulated up to tmax- A single run is noisy and one 
has to repeat many times, order 10^. The much more common choice is to start from a 
reasonable nonequilibrium configuration and to equilibrate before measuring correlations. 
Usually one then simulates very long trajectories, up to times of order 2^®, over which the 
time lag ga{j,t + 'r) 5 'o'( 0 ,t) is sampled, r < tmax- In addition one averages over a small 
number of runs, of order 10^. The total number of samples is roughly the same in both 
methods. The random number generator method produces the thermal average with a 
higher reliability. 

The sampled Saa'U, t) can be Fourier transformed in the spatial variable and/or in the 
time variable. One can also transform to normal coordinates. These are linear operations 
which can be done for each sample or only after averaging. Should one keep the full 
resolution or only some data points? Of course it depends somewhat on the goals. In 
[58, 59], the full 3x3 matrix is sampled and subsequently transformed via the theoretically 
computed R matrix to obtain This approach allows to test diagonality. Because 

the peak structure is most easily seen in the space-like j coordinate, maximal resolution 
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Figure 3: Scaling plot of heat and sound peak for a FPU chain with N = 8192, potential 
parameters a = 2, (3 = 1, pressure -P = 1, and temperature = 0.5. 

for j is retained and only three times {t = 250, 500,1000) are recorded for the purpose 
of making a scaling plot. In [63], the lowest Fourier modes are measured as a function 
of t. In [57] the lowest Fourier modes are plotted as a function of the frequency u. A 
separate issue are the much simulated total current correlations. The total currents are 
sampled directly, in the most complete version momentum, energy, and cross correlations, 
and then plotted as a function of f or cj. 

We reproduce only a few hgures. Many more details can be found in the original 
papers. In Fig. 3 we display the data for the FPU chain [58] with lA’pu(a:) = \x^ + |x 3 + 
/d = 2, P = 1, 77 = 8192, and fmax = 2700. The sound speed is c = 1.45. Note that 
the sound peak is somewhat distorted, not symmetric relative to ct, but has a rapid fall-off 
away from the sound cone. As only £t parameter one uses Ah, resp. Ag. The optimal £t 
at the longest available time is denoted by A]]™'’, resp. Ag™P, standing for empirical value. 
In most cases there is also a theoretical value based on decoupling and/or mode-coupling, 
which is indicated in square brackets. For the FPU simulations the results are for the 
heat peak A|]™^ = 13.8 [1.97], and for the sound peak A|™p = 2.05 [0.68]. From the visual 
appearance, one might have guessed the theoretical values to be just the other way round. 
So maybe the system tries to generate the optimal Levy peak at non-relaxed sound peaks. 

In Fig. 4 we reproduce the plot of heat and sound peak for the hard-collision models 
with shoulder, hard-point, and square well potential, in the latter two cases with alternat¬ 
ing masses [59]. To better judge the quality of the fit we provide a logarithmic plot. In all 
cases N = 4096 with = 1024. The following parameters have been chosen, shoulder: 
Co = 1, P* = 1.2, /3 = 2, c = 1.74, hard-point: mi/rriQ = 3, P = 2, /d = 0.5, c = 1.73, 
square well: mi/rriQ = 3, a = 1, P = 0, /3 = 2, c = 1.73. The £t to the predicted scaling 
function has an error less than 5%. For the hard-collision models the simulation results 
are, shoulder: Af“P = 1.62 [1.71], = 1.44 [1.04], hard-point: A|“p = 1.42 [2.00], 

Aemp _ [0.95], square well: A^™'’ = 0.95 [1.04]. Recall that the square well potential 
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Figure 4: Heat and sound peak for shoulder, hard-point, and square well potential in 
logarithmic scale. 


at P = 0 is in a distinct dynamical universality class, hence the different scaling expo¬ 
nents. For this model As is related to a diffusion constant, which can be obtained only 
numerically. 

The MD simulation [57] is at low temperatures and considers the positional correla¬ 
tions, which in Fourier space differ from the stretch correlations only by a fc-dependent 
prefactor. Thus sampled is the correlator 

(1 — cos(27rfc))“^ y dte^'^^Sii{k,t) = {1 — cos(27ik))~^Sii{k,u;). (5.4) 

For low temperatures the area under the central peak is a factor 10 smaller than the one 
under the sound peak. Hence only the sound peak is explored. Its asymptotic scaling 
form is 

Sn{k, t) = cos{2Trikct)CnfKPzik{Xh\t\f^^). (5.5) 

Considering only the right moving sound peak by setting (Umax = 27ikc, 

Su{k,u + u^^,) = j 

= J (5.6) 

Thus defining 

hKPz(i^) = J dte^‘^^fKPz{\k ^^)! (5-7) 
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one arrives at 


(5.8) 


o; + Wmax) — ^/iKPz(l^/All|fcp/^) . 

If one normalizes the maximum of S to 1, then hKPz(i^) is replaced by hKPz(i^)/hKPz(0), 
which amounts to setting the prefactor in (5.8) equal to 1. In Fig. 5 the spectrum at two 
different choices of the asymmetry parameter is displayed. 

In the simulation the potential is chosen as V{x) = + \x'^, where the 

asymmetry varies from 0 to 2. Increasing a, the pressure P increases from 0 to 0.2 and 
the inverse temperature from 9.55 to 9.75. The sound speed c ~ 1.1. In frequency space 
the peak moves linearly in k at around u = 0.01. For the shape function one uses S{k,u}) 
for k = 1,2,4, 8,16 to generate a scaling plot. Over the whole range of a’s the £t with 
the scaling function (5.7) is fairly convincing. The optimal fit parameter starts from 
yemp ^ q_q 2 [Q.Od] at « = 0.2 to A®”P = 0.08 [0.37] at a = 1.6 and to = 0.07 [0.53] 
at a = 2.0. 

KPZ scaling - FPU a(3 


e=0.5 a=0.1 P=1 



Figure 5: Scaling plot of the sound peak in {k = q,u}) variables of a FPU chain with 
N = 1024, potential parameters a = 0.1, /3 = 1, pressure P = —0.04, and temperature 
= 0.105. The speed of sound is c = 1.11. 

Given the diversity of models, parameters, and numerical schemes only tentative con¬ 
clusions can be drawn. 

(i) The separation into three peaks is a fast process. In normal mode representation 
the off-diagonal matrix elements are indeed small. There are no correlations beyond the 
sound cone. 

(ii) The distinction between G\i ^ 0 and G\i = 0 is seen very convincingly. 
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(in) The central peak is on the scale and adjusts well to the predicted Levy dis¬ 
tribution with one caveat. As seen from the center, the shape function fairly rapidly 
switches into the power law decay. However at the location of the sound peaks there are 
some wiggles and beyond only small amplitude noise is observed. In this sense, the Levy 
distribution gets uncovered as time progresses. The Levy distribution is a fairly direct 
consequence of Gqq = 0 and thus one of the strongest supports for the theory. 

(iv) The sound peaks are mirror images of each other. One plots them on the scale, 
but then there is still a slow change in time. For the hard-collision models the shape is 
almost perfect, but the As parameter is dropping in time. For the FPU models the peak 
is distorted and still away from the symmetric shape predicted by the theory. To the 
outside of the sound cone there is the rapid fall-off in accordance with the KPZ scaling 
function. But towards the heat peak there is slow decay. The sound peak is tilted away 
from the central peak. Apparently there is still a strong interaction between the peaks. 

(v) \n. the majority of the simulations the peaks vary slowly on the scale for sound, 

resp. for the heat peak. Thus it becomes meaningful to use as a £t the theoretical 
scaling function with Ag, resp. Ah, as only free parameter. This optimal choice has been 
denoted above as empirical value The error between the measured and theoretical 

shape function is less than 5%. However one observes that Ag™P and A^™’’ are still changing 
in time signaling that the simulation has not yet reached the truly asymptotic regime. In 
some simulations A^^^ drops monotonically in time and differs not too strongly from Ag/h- 
One is then willing to believe that for even longer simulation times the asymptotic value 
is reached. But in other simulations there is a much stronger discrepancy, which asks for 
more explanations. 

Mode-coupling is not specihc to anharmonic chains. In principle any one-dimensional 
system with conserved helds can be handled by the same scheme. This offers the possibility 
to test the theory through other models, possibly finding systems with less strong hnite 
time effects. One obvious choice are stochastic lattice gases with several type of particles 
like several lane TASEP [64, 43] and the AHR model [65, 66]. For them the couplings can 
be more easily adjusted than for anharmonic chains, which offers the possibility to test the 
dynamical phase diagram. Also anharmonic chains with a stochastic collision mechanism, 
respecting the conservation laws, have been studied in considerable detail [67, 68, 42]. 

6 Total current correlations 

The total current is a fluctuation observable, in contrast to Saa'{j,t) which refers to the 
average of the product of two local observables. Thus we need some additional consider¬ 
ations to establish the link to nonlinear fluctuating hydrodynamics. For a ring of size N, 
A]y = [1,..., 77], the total currents are dehned by 

1 ^ ^ 

( 6 . 1 ) 
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and the total current covariance reads 


N 

rAjv,a«'(^) (t7tot,Ajv,o(^)) •-^ot,Ajv,o' (0)) P,I3An P,I3,An ■ (6-2) 

t=l 

The cumulant {■; ■) means that the static average is subtracted and system size is indicated 
explicitly. In the limit N ^ oo 

^aa'it) = ^(j7Q,(j, t); J7a'(0, 0))p,/3 . (6.3) 

For hxed t the integrand decays exponentially in j, but with a correlation length increasing 
in time. 


r22(t) 



(a) momentum current correlations 


r22(t)-OT chy6 



r33(t) 



(c) energy current correlations 





Figure 6: Total momentum and energy correlations for a hard-point particles with alter¬ 
nating masses 


Before we emphasized the transformation to normal coordinates. For the currents we 
stick to the physical helds. The current j7tot,Aiv,i(^) itself conserved. Thus only the (2, 3) 
block has a time variation. Using stationarity and time-reversal, the diagonal elements 
are even, Taa(t) = Taa{—t), while the off-diagonal elements are odd and satisfy 

r23(^) = —r23(—= r32(—t) = —F32(f) . (6.4) 
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Mode-coupling predicts that this matrix element vanishes. In our simulations we observe 
an exponential decay with a decay time of order 20 to 30. Thus the correlations of real 
interest are T 22 {t) and r 33 (t). For chains the latter is the most frequently simulated 
equilibrium time correlation. The momentum current correlations have been measured in 
[53, 69] and the momentum-energy cross correlations only recently in [70]. 

rA]v,aa'(^) is a fluctuation observable, for which the equivalence of ensembles does not 
hold. But the time-independent difference between microcanonical and canonical average 
can be computed explicitly. For the microcanonical ensemble, limi^oo Taj^, micro(t) = 0. 
On the other hand, there is no reason for to vanish asymptotically if the canonical 

average, as in (6.2), is used. In fact, for inhnite volume, 

lim F 22 (t) =, lim F 33 (f) =. (6.5) 

t^OO t^OO 

The asymptotic values in (6.5) are called Drude weight, which has received a lot of 
attention in the context of current correlations for integrable quantum chains. A non¬ 
zero Drude weight indicates that the correlator of the corresponding conserved held has a 
ballistically moving component, but it cannot resolve the structure of this component. For 
non-integrable anharmonic chains, the ballistic pieces are just the sharply concentrated 
sound peaks, while in the integrable case one expects to have a broad spectrum which 
expands ballistically. 

The link to mode-coupling is achieved through the general observation that the current 
correlations are proportional to the memory kernel. We use its diagonal approximation 
and insert the asymptotic form of fa- The memory kernel is in normal mode representa¬ 
tion. Thus we still have to transform back to the physical helds through the R matrix. 
The computation can be found in [70] with the result 

r^ 2 (t) = F 22 (t)-rV 

^ |(Aht)“^/®('i/’ 0 ,^f>o)^y dx/Levy,5/3(3;)^ + (Ast)“^/^(^/’l,i 7 >l)^y' dx fKPz{xY 

= F 33 (t) - ~ c^/3"^(Ast)"^/^ j dx /kpz(x)^ , (6 .6) 

where {fja} are the eigenvectors of A, = 0, = cfji, see [29]. If (4.33) is satished, 

e.g. an even potential at P = 0, then the sound peak is diffusive, see Eq. (4.35), and the 
central peak is Levy 3/2. Furthermore ('/05 = 0 = (-/i, Hence T22{t)—(3~^c^ 

is expected to decay integrably, while 

^ c2/3-2(87rAf)-'/2 _ 

The energy current correlation is predicted to decay as which has been reported 
in MD simulations already more than 15 years ago [51, 52]. The true mechanism behind 
the decay is actually somewhat subtle. From the conservation law it follows that the 
second moment of the heat peak is related to the second time derivate of the current 
correlation. Using the asymptotic form (4.27), including the cut-off at the sound peak, 
one arrives at 

[ dxx‘^fo{xR) ~ |^(Ah)®/^c^''^(Aht)“^/^ 

J —ct 


dR 
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( 6 . 8 ) 


This argument overlooks that the scaling form is for the normal mode representation, while 
r 33 (t) refers to the physical energy current. The complete computation leads however to 
the same power law except for a different prefactor [70]. 


r22(t) 



(a) momentum current correlations 


T22{t)-C^IIS 




Figure 7: Total momentum and energy current correlations for a hard collision model 
with shoulder potential. 

The momentum current correlation should decay as which is a recent hnding. 

However, its prefactor could vanish, in principle. Mode-coupling with the currently avail¬ 
able precision would not provide an answer, then. In contrast the prefactor for r 33 (t) is 
strictly positive. 

Our simulation results [70] are shown in Fig. 6,7. The parameters are as before, 
N = 4096, fmax = 1024, shoulder: cq = 1, T* = 1.2, (5 = 2, c = 1.74, and hard-point: 
mi/mo = 3, P = 2, (3 = 0.5, c = 1.73. The red lines indicate the predictions based on 
mode-coupling. It is interesting to note that for the shoulder potential the evidence for 
a decay is not so overwhelming as one might have anticipated and, by looking at 
a different time window, one could as well £t to a slightly different exponent. On the 
other hand the hard-point potential with alternating masses shows a very clean power 
law decay. Through MD with shoulder potential the predicted decay of T 22 {t) is well 
conhrmed. However, for the hard-point potential it so happens that both prefactors, 
{'ipOyH'^'ipo) and {'ipi, vanish. Numerically we estimate a decay as Again this 

is a strong qualitative support of nonlinear fluctuating hydrodynamics. One might have 
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thought that all hard-collision potentials have the same asymptotic power law for the 
momentum current correlation. This expectation is born out under the proviso that the 
respective prefactors do not vanish. Matrix elements, as {ijjo, , must come from a 

microscopic computation and cannot be deduced by a mere inspection of the potential. 

An additional conhrmation of (6.5), (6.7) has been accomplished recently [71]. For 
the potential V (x) = /^ = 1 and at P = 0.59 with a = —2, resp. at P = —0.5, 

a = 1.89, up to very small errors the signature of the G matrices is identical to an even 
potential at P = 0. In the MD simulation the energy current is found to decay as 
and the momentum current seems to be integrable. Keeping all parameters hxed and 
shifting slightly to a = —2.7, resp. a = 2, the decay as stated in (6.6) is restored. 

7 Other ID Hamiltonian systems 

For the anharmonic chains studied so far, the potential depends only on qj+i—qj and hence 
remains without change under spatial translations. Physically this property is obvious 
and seems hard to avoid. However there could be a substrate potential which forces 
the particles preferentially to particular locations. One could consider a two-component 
system, which then has an acoustic and an optical mode. The latter would be comparable 
to a one-component system with an on-site potential. Such considerations lead to the more 
general class of Hamiltonians 

N N-1 

»»= E (h?++ E - ft) P-1) 

i=i i=i 

with some conhning on-site potential I4s. The only conserved held is the energy. There 
is a unique equilibrium measure. The Euler currents vanish. From the perspective of 
huctuating hydrodynamics all evidence points towards diffusive energy transport. For the 
case of a quadratic V and for = Pppu very detailed MD simulation conhrm diffusive 
transport [72]. 

More interesting are models with two conservation laws. We discuss separately coupled 
rotators, which can be thought of as a classical limit of a quantum Heisenberg chain, and 
the discrete nonlinear Schrodinger equation on a lattice, which is the classical held theory 
for lattice bosons. 

Coupled rotators. The Hamiltonian of the rotator chain reads 

N 

= E (M + - Vi)) (7.2) 

i=i 

with periodic boundary conditions, (pw+i = At hrst glance we have only rebaptized 
qj as ifj. But the (pjS are angles and the p^-’s angular momenta. Hence the phase space is 
{S^ X M)^ with denoting the unit circle. The standard choice for V is K('d) = — cos'd, 
but in our context any 27r-periodic potential is admitted. The equations of motion are 
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Obviously angular momentum is locally conserved with the angular momentum current 


Mi) = - iPj-i). 


(7.4) 


As local energy we dehne Cj = Then ej is locally conserved, 

= pj+iV'{ipj+i - pj) - PjV'ipj - , 

from which one reads off the energy current 

J2U) = -PjV'i^j - Pj-i). 


since 


(7.6) 


(7.6) 


For the angles, in analogy to the stretch, one dehnes the phase difference fj = 

where 0 is 27r-periodic and 0(x) = x for |x| < tt. Because of the jump discontinuity the 

stretch is not conserved. A rotator chain has only two conserved helds. 

To apply nonlinear fluctuating hydrodynamics one has to compute the Euler currents 
in local equilibrium. Since there are two conserved helds the canonical equilibrium state 
reads 


1 ^ 

— JJexp [ - P{^{pj - uf + V{pj+i - pj))]dpjdpj (7.7) 

with u the average angular momentum. Now 

{Jlij))N = —{y'iPj — Pj-l))N , {J2{.3))n = —u{V\(pj — , (7.8) 

average with respect to the canonical ensemble (7.7). We claim that 

lim {V\pj - Pj-i))N = 0 . (7.9) 

A'^OO 

For this purpose we expand in Fourier series as 

g-vW = , /(i9) = ^ . (7.10) 

mGZ mGZ 

Then, working out all Kronecker deltas from the integration over the (p/s, one arrives at 
(/((y9j+i - V?i))iv = ( a(m)^) ^a(m)^(^a(m)-i^/(£-m)a(£)) . (7.11) 


mez 


rnGl, 


fez 


Since a(0) > |a(m)| for all m 7 ^ 0, 

lim {f{(pj+i - (pj))N = a( 0 )”^ 


(7.12) 


For /('d) = W('d), the latter integral vanishes because of periodic boundary conditions in 
We conclude that both currents vanish on average. 
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As before we consider the infinite lattice with thermal expectation and form the 
eqnilibrinm time correlations as 

Sac'U.t) = - (^a(0))„,/3(^o(0))„,/3, (7.13) 

g{j,t) = [pj{t),ej{t)y pj is odd and Cj is even nnder time reversal. Hence in the Green- 
Knbo formnla the cross term vanishes. Thns, for large j, t, flnctnating hydrodynamics 
predicts 

Saa’ij,^) = , (7.14) 

where /g is the nnit Ganssian. Da is the diffnsion coefficient of mode a. Of conrse, it can 
be written as a time-integral over the corresponding total cnrrent-cnrrent correlation, bnt 
its precise valne has to be determined nnmerically. This has been done for the standard 
choice V{x) = — cosx, to which we specialize now. For /3 = 1, with a lattice size N = 500, 
the diffnsive peaks are well established at f = 2000 [73, 74]. Energy diffnsion has been 
conhrmed mnch earlier [75, 76]. 

At low temperatnres one hnds a different, perhaps more interesting scenario. At zero 
temperatnre, there is the one-parameter family of gronnd states with ipj = p, pj = 0. 
When heating np, nnder the canonical eqnilibrinm measnre, the phase pj jnmps to Pj+i 
with a jump size 0{l/^/^). Next we have to nnderstand how the conservation of rj held is 
broken. In a pictorial langnage, the event that \pj^i{t) — Pj(t)\ = tt is called an nmklapp 
for phase difference fj or an nmklapp process to emphasize its dynamical character. At 
low temperatnres a jnmp of size tt has a small probability of order with AH = 2 

the height of the potential barrier. Hence fj is locally conserved np to nmklapp processes 
occnrring at a very small freqnency only. This can be measnred more qnantitatively by 
considering the average 

ruk(t) = ((G(t)7o(0))n,/3 - • (7.15) 

iez 

At /3 = 1, ruk(t) decays exponentially dne to nmklapp. Bnt at /? = 5 the decay rate is 
already very mnch snppressed, see [74]. 

In the low temperatnre regime it is tempting to nse an approximation, where the 
potential V{x) = — cosx is Taylor expanded at the minimnm x = 0. Bnt snch procednre 
wonld nnderestimate the regime of low temperatnres, as can be seen from the example of 
a potential, still with AH = 2, bnt several shallow minima. The proper small parameter 
is snch that /dAH > 1. To arrive at an optimal low temperatnre hamiltonian, we 
hrst parametrize the angles (pi,...,(pAr throngh Vj = Pj+i — pj with rj G [—7r,7r]. To 
distingnish, we denote the angles in this particnlar parametrization by (pj. The dynamics 
governed by TTcr corresponds to periodic bonndary conditions at Vj = ivr. For a low 
temperatnre description we impose instead specnlar reflection, i.e., if Vj = ivr, then 
Pj, pj+i are scattered to p'j = Pj+i, p)+i = Pj. By hat all nmklapp processes are now 
snppressed, while between two nmklapp events the GR dynamics and the low temperatnre 
dynamics are identical. The corresponding hamiltonian reads 

N 

HcRjt = (p? + (7.16) 

i=i 
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with 

y(a;) = —cosx for |a:|<7r, y(a:) = oo for |a;|>7r, (7-17) 

periodic boundary conditions 0Ar+i = 0i being understood. The pair {(j)j,Pj) are canon¬ 
ically conjugate variables. Note that as weights exp[—/^ifcn] = exp[—Thus all 
equilibrium properties of the coupled rotators remain untouched. 

The hamiltonian hfcR,it is a variant of the hard collision model with square well poten¬ 
tial as discussed before, see [59]. The dynamics governed by hfcR.it has three conserved 
fields, the stretch rj = 0^+1 — cjyj, the momentum and the energy Cj = -|- V{rj). 

Because of 0i = (j)N+i, one has — 0- The model is in the dynamical phase 

characterized by an even potential at zero pressure. 

We claim that, for 0A1/ > 1, the CR equilibrium time correlations are well approx¬ 
imated by those of hfcR.it, provided the time of comparison is not too long. The latter 
correlations can be obtained within the framework of nonlinear fluctuating hydrodynam¬ 
ics. Thereby one arrives at fairly explicit dynamical predictions for the low temperature 
regime of the CR model. 

One physically interesting information concerns the Landau-Placzek ratio at low tem¬ 
peratures. We use (4.29) and expand R~^ in 1/0. To lowest order it suffices to use the 
harmonic approximation, V(x) = ^ax^, a = 1 for the cosine potential. Using the formulas 
in Appendix A of [29] one arrives at the following value for the Landau-Placzek ratios, 

r-r : (2a0)“^(l, 0,1) , p-p : (20)“^(1,0,1), e-e : (20)“^(a£^, 0“^, af'^). (7.18) 

The correlations are small, order 0“^. For the stretch correlations there is no central 
peak, to this order, and for the energy correlations the central peak is down by a factor 
0 “^ relative to the sound peaks. 

To have a unihed picture we add fj to the list of helds of physical interest. At high 
temperatures rj is not conserved and one has diffusive spreading of the conserved helds. 
At low temperatures fj is conserved up to small errors and the conventional three-peak 
structure, including universal shape functions, results. For extremely long times umklapp 
processes will happen and one expects that they force a cross over to the Gaussian scaling 
(7.14). The precise dynamical structure of such cross over still needs to be investigated. 

Nonlinear Schrddinger equation on a lattice. A further example with a dynam¬ 
ically distinct low temperature phase is the nonlinear Schodinger equation on the one¬ 
dimensional lattice. In this case the lattice held is ipj G C, for which real and imaginary 
part are the canonically conjugate helds. The Hamiltonian reads 

N 

^ = ('^•19) 

i=i 

with periodic boundary conditions and coupling g > 0. The sign of the hopping term plays 
no role, since it can be switched through the gauge transformation 'ijjj The chain 

is non-integrable and the locally conserved helds are the number density pj = 10^0 and 
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the local energy e^- = — V'iP + Hence the canonical equilibrium state is 

given by 

N N 

JJ ^ N = iV'iT , (7.20) 

i=i i=i 


with the chemical potential /i. We assume (3 > 0. But also negative temperature states, in 
the microcanonical ensemble, have been studied [77, 78]. Then the dynamics is dominated 
by a coarsening process mediated through breathers. In equilibrium, the "^-held has 
high spikes at random locations embedded in a low noise background, which is very 
different from the positive temperature states considered here. For them the density 
and energy currents are symbolically of the form i{z — z*), hence their thermal average 
vanishes. Both fields are expected to have diffusive transport. In fact, this is confirmed 
by MD simulations [79]. They also show Gaussian cross-correlations, which is possible 
since density and energy are both even under time reversal. In the previous studies [80] 
transport coefficients have been measured in the steady state set-up. 

To understand the low temperature phase, it is convenient to transform to the new 
canonical pairs Pj,(pj through 

V., = v^e‘« . (7.21) 


In these variables the Hamiltonian becomes 


N 

H = '^{- pj cos(<^j+i - ifj) + pj + \gp]) . (7.22) 

t=i 


The equations of motion read then 


dtp, = -dp^H , dtp, = d^^H . (7.23) 

Pj takes values on the circle and pj > 0. From the continuity of ipjit) when moving 
through the origin, one concludes that at pj(t) = 0 the phase jumps from Pj(t) to pj(t)+7i. 

One recognizes the similarity to the coupled rotators (7.2). But now the equilibrium 
measure carries a nearest neighbor coupling. For /x > 0, in the limit /? —)■ cx) the canonical 
measure converges to the one-parameter family of ground states with pj = p = p/g, 
Pj = p with p uniformly distributed on S^. At low temperatures the field of phase dif¬ 
ferences Tj is approximately conserved. The low temperature hamiltonian is constructed 
in such a way that the equilibrium ensemble remains unchanged while all umklapp pro¬ 
cesses are suppressed. To achieve our goal we follow verbatim the CR blueprint. The 
phases are parametrized such that Pj+i — Pj lies in the interval [—vr, tt] and this particular 
parametrization denoted by (pj. Umklapp is a point at the boundary of this interval. Now 
{pj, Pj) are a pair of canonically conjugate variables, only pj > 0 instead of pj G M. Thus 
the proper low temperature hamiltonian reads 

N-l 

Hit = (vPi+i pj U{pj+i - Pj) + u(Pi)) , (7.24) 

j=0 
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where 


U{x) = —cos{x) for |a:|<7r, U{x) = oo for |x|>7r, (7.25) 

and 

V{x)=x+^gx‘^ for x>0, V{x) = oo for x<0. (7.26) 

The low temperature Hamiltonian has a nearest neighbor coupling, which complicates 
the scheme through which the G matrices are determined [79]. Progress is achieved 
through the miraculous identity j = (/x, P, fiP) for the Euler currents. The G coefficients 
are evaluated at P = 0. As for a generic anharmonic chain, one hnds that Gqq = 0 and 
GJi 7 ^ 0. Thus the heat peak is predicted to be Levy 5/3 and the sound peaks to be 
KPZ. The sound peak for the density-density correlation was hrst observed in [81] using 
k,uj space, see also [82]. In [79] we use normal mode representation, as explained in this 
article. The sound peaks £t nicely with KPZ, but the normalized heat peak is very broad 
and noisy, still with a shape not unlikely Levy 5/3. 
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